% This script rescales and plots impulse response functions for 
% contemporaneous or news TFP shocks.

% Christoph Gortz


close all
clear all

%% Initialisation and load data

load Output_BayesianIRF_new2
%load BayesianIRF_TFPobs     % this is output from estimation with TFP as observable

shockPlot = 'z';              % Shock process to be plotted
exo_shockPlot = 'ez8';        % Exogenous shock to be plotted
shockOther = 'v';             % Other Shock process with trend
exo_shockOther = 'ev';        % Other Exogenous shock with trend

% Variables names from dynare to plot
vars = [
    'y         '; 
    'c         '; 
    'r         '; 
    'sc        '; 
    'si        '; 
    'spreadobsc';  
    'inve      '; 
    'invec     '; 
    'invei     '; 
    'nc        '; 
    'ni        '; 
    'spreadobsi';  
    'lab       '; 
    'labc      '; 
    'labi      '; 
    'lrc       '; 
    'lri       '; 
    'z         '; 
    'pinfc     '; 
    'pinfi     '; 
    'w         '; 
    'qc        '; 
    'qi        '; 
    'pi        '; ];
% Names of subplots
IRFtitle = [ 'Output                   ';
            'Consumption              ';
            'Nom. Interest Rate       ';
            'C-Sector Capital         ';
            'I-Sector Capital         ';
            'C-Sector Spread          ';
            'Total Investment         ';
            'C-Sector Investment      ';
            'I-Sector Investment      ';
            'C-Sector Bank Equity     ';
            'I-Sector Bank Equity     ';
            'I-Sector Spread          ';
            'Total Hours              ';
            'C-Sector Hours           ';
            'I-Sector Hours           ';
            'C-Sector Leverage Ratio  ';
            'I-Sector Leverage Ratio  ';
            'C-Sector TFP             ';
            'C-Sector Inflation       ';
            'I-Sector Inflation       ';
            'Real Wage                ';
            'C-Sector Price of Capital';
            'I-Sector Price of Capital';
            'Relative Price of Inve.  ';
            ];       
nVars = size(vars,1);

figDimensions = [4, 6];     % Dimensions of figure - needs to be consistent
                            % with nVars
irfLength = 40;             % Length of generated IRFs 

FileName = ['BayesianIRF_', exo_shockPlot];


%% LOAD INFO AND DETREND

% Stores lower percentile, mean and upper percentile of Bayesian IRFs for
% each variable specified in vars
for iVars = 1:nVars
    var = strcat(vars(iVars,:));
    eval([ 'percentiles =  oo_.PosteriorIRF.dsge.Distribution.' var '_' exo_shockPlot ';']);
    %eval([ 'percentiles =  oo_.PosteriorIRF.dsge.deciles.' var '_' exo_shockPlot ';']);
    lowerPercentile = percentiles(:,1);
    upperPercentile = percentiles(:,9);
    eval([ 'mean =  oo_.PosteriorIRF.dsge.Mean.' var '_' exo_shockPlot ';']);
    eval([ 'BayesIRF_' var ' = [lowerPercentile , mean, upperPercentile] ;']);
end


% Load shock's irf data
eval([ shockPlot 'Mean = oo_.PosteriorIRF.dsge.Mean.' shockPlot '_' exo_shockPlot ';' ])
eval([ shockOther 'Mean = oo_.PosteriorIRF.dsge.Mean.' shockOther '_' exo_shockPlot ';' ])


% Transform series of shock z to isolate zd:
vd =[]; vp = 0;
zd =[]; zp = 0;

for m = 1:numel(zMean)
    % Change of v and z in percent
    vp = vMean(m,1)+vp;
    vd = [vd; vp];
    zp = zMean(m,1)+zp;
    zd = [zd; zp];
end    
vd = [vd, vd, vd];
zd = [zd, zd, zd];




%% Descale stationarised variables: 
% According to normalisation given fore each variable, e.g. K_t = K_t^~*Z_t
% Variables need to be called like in variable 'vars' with addition
% '_detrend' for plotting
y_detrend = BayesIRF_y + (alfac/(1-alfai))*vd + zd;
c_detrend = BayesIRF_c + (alfac/(1-alfai))*vd + zd;
spreadobsc_detrend = BayesIRF_spreadobsc;
spreadobsi_detrend = BayesIRF_spreadobsi;
inve_detrend = BayesIRF_inve + (1/(1-alfai))*vd;
invec_detrend = BayesIRF_invec + (1/(1-alfai))*vd;
invei_detrend = BayesIRF_invei + (1/(1-alfai))*vd;
qc_detrend = BayesIRF_qc + zd - ((1-alfac)/(1-alfai))*vd;
qi_detrend = BayesIRF_qi + zd - ((1-alfac)/(1-alfai))*vd;
lab_detrend = BayesIRF_lab;
labc_detrend = BayesIRF_labc;
labi_detrend = BayesIRF_labi;
pinfc_detrend = BayesIRF_pinfc;
pinfi_detrend = BayesIRF_pinfi;
% zcapcd = Variables{12,1};
% zcapid = Variables{13,1};
% rkcd = Variables{14,1} + zd - ((1-alfac)/(1-alfai))*vd; % Is this correct?
% rkid = Variables{15,1} + zd - ((1-alfac)/(1-alfai))*vd;
% mccd = Variables{16,1};                                 % Is this correct?
% mcid = Variables{17,1};
w_detrend = BayesIRF_w + zd + alfac/(1-alfai)*vd;
lrc_detrend = BayesIRF_lrc;
lri_detrend = BayesIRF_lri;
nc_detrend = BayesIRF_nc + (alfac/(1-alfai))*vd + zd;
ni_detrend = BayesIRF_ni + (alfac/(1-alfai))*vd + zd;
r_detrend = BayesIRF_r;
sc_detrend = BayesIRF_sc + (1/(1-alfai))*vd;
si_detrend = BayesIRF_si + (1/(1-alfai))*vd;
%z_detrend = BayesIRF_z;
z_detrend = BayesIRF_z;
pi_detrend = BayesIRF_pi + zd + (alfac-1)/(1-alfai)*vd;

%% Plot IRFs and save
figure1 = figure('Color',[1 1 1]); % Generates white background around plot
for iVars = 1:nVars
        plotVar = strcat(vars(iVars,:));
        irfTitle = strcat(IRFtitle(iVars,:));
        % Subplots
        subplot(figDimensions(1),figDimensions(2),iVars);
        eval([ 'h = area([' plotVar '_detrend(:,1), (' plotVar '_detrend(:,3)-' plotVar '_detrend(:,1)) ]);' ]);
        set(h(2),'FaceColor',[.9 .9 .9]); % Shades the area between curves
        set(h(1),'FaceColor',[1 1 1]);    % White background in plot
        eval([ 'set(h,''BaseValue'',min(' plotVar '_detrend(:,1)));' ]); % Excludes line at zero
        hold on
        eval([ 'plot(' plotVar '_detrend(:,2),''-k'',''linewidth'',3);' ]);
        axis tight
        hold off
        eval([ 'title( irfTitle );' ]);
end


